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Cygnus X-1 contains a 21-solar mass black hole — 
implications for massive star winds 
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The evolution of massive stars is influenced by the mass lost to stellar winds 
over their lifetimes. These winds limit the masses of the stellar remnants (such 
as black holes) that the stars ultimately produce. We use radio astrometry to 
refine the distance to the black hole X-ray binary Cygnus X-1, which we find 
to be 2.22*015 kiloparsecs. When combined with previous optical data, this 
implies a black hole mass of 21.2 + 2.2 solar masses, higher than previous mea- 
surements. The formation of such a high-mass black hole in a high-metallicity 
system constrains wind mass loss from massive stars. 


Gravitational wave detections of black hole merger events have revealed a population of 
black holes with masses ranging from 7 to 50 solar masses (Mo) (1). Black holes that interact 
with a companion star are visible to electromagnetic observations as an X-ray binary. Radial 
velocity measurements of these companion stars have shown that black holes in X-ray binaries 
all have masses below 20M. (2). The highest measured black hole mass in an X-ray binary is 
15.65 + 1.45 M; for the extragalactic system M33 X-7 (5). 

The mass of a black hole is initially set by the properties of its progenitor star, augmented by 
accretion or mergers over its lifetime. The relevant properties of the progenitor include its initial 
mass and abundance of heavy elements (referred to as its metallicity), the mass lost in stellar 
winds over its lifetime, and the evolutionary pathway that it followed, which can be strongly 
influenced by a binary companion. Mass measurements for massive stellar-mass black holes 
constrain stellar and binary evolution models (4) and predictions of the expected black hole 
merger rates. 

The X-ray binary Cygnus X-1 (V1357 Cyg; coordinates in Table S1) contains a black hole in 
a 5.6-day orbit with a more massive supergiant donor star, of spectral type O. Previous estimates 
of its component masses were based on a parallax measurement, whereby the apparent annual 
angular shift in the source position relative to more distant objects was measured using radio 
very long baseline interferometry. Via trigonometry, this gave a distance of 1.8670 1; kpe (5). 
When combined with optical radial velocity measurements of the system, these data yielded a 
black hole mass of 14.843- 1.0M 5 (6). However, the derived system parameters are then inconsis- 
tent with the expected mass-luminosity relation for the donor, if it is a hydrogen-burning main 
sequence star (7). The optical parallax measurement of 0.42 + 0.03 milliarcseconds (mas) using 
the Gaia space telescope (6), after correction for the known zero-point offset of ~ 0.05 mas 
(with estimates ranging from 0.03—0.08 mas; (9)), becomes 0.47 + 0.04 mas. This is inconsis- 
tent with the radio value of 0.54 +0.03 mas (5), and is unlikely to be due to orbital displacement 
of the donor star because the Gaia value is the average over 119 orbital periods. 

Between May 29 and June 3, 2016, we performed six observations (one per day) of Cygnus X- 
1 with the Very Long Baseline Array (VLBA) at 8.4 GHz, sampling one full orbital period. To 
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reduce systematic uncertainties in our position measurements, we phase referenced the data to 
a nearby calibrator source 0.4? from Cygnus X-1 (70). These data resolve the orbital motion 
of the black hole as projected onto the downstream surface from which the jet emission can es- 
cape, known as the photosphere. We find the orbit is clockwise on the plane of the sky (Fig. 1), 
in agreement with previous observations (5). 

Combining the orbital phase coverage of our VLBA data with the archival observations (5), 
we simultaneously fitted the entire data set (covering a 7.4-year baseline) with a full astrometric 
solution incorporating linear motion across the sky (proper motion), parallax and orbital motion. 
A fit to both right ascension and declination co-ordinates showed an orbital phase dependence in 
the direction of the residuals. We attribute this to the effect of free-free absorption in the stellar 
wind, which is known to modulate the radio emission of Cygnus X-1 on the orbital period (//). 
Electrons in the ionised wind of the O star can absorb radio photons in the presence of atomic 
nuclei, preventing the radio emission from the inner parts of the jet from reaching us. This free- 
free absorption is reduced as the stellar wind density decreases on moving away from the O star, 
allowing radiation to escape from further downstream. As the black hole (from which the jet is 
launched) moves around its orbit, the varying path length through the stellar wind imprints an 
orbital periodicity on the apparent radio position along the jet axis. When the black hole is on 
the far side of the donor star, the path length and hence absorption are maximised, pushing the 
radio photosphere downstream along the jet axis (70). 

To negate the impact of the stellar wind absorption, we therefore repeated our astrometric 
model fitting in one dimension only, perpendicular to the known jet axis. This removed the 
orbital phase dependence of the fit residuals perpendicular to the jet axis (Fig. 2). Our measured 
semi-major axis of the black hole orbit is 58 + 20 microarcseconds (uas), and our revised 
parallax measurement is 0.46 + 0.04 mas, consistent with the optical value from Gaia after 
correction for the zeropoint. After converting our measured parallax to a distance using an 
exponentially-decreasing space density prior (/2), we find a distance of 2.221018 kiloparsecs 
(kpc). 

This revised distance affects the system parameters derived from the optical modelling (6). 
We reanalysed the optical light curve (/3) and radial velocity curves (/4), adopting our revised 
distance and additional constraints on the effective temperature, surface gravity and helium 
abundance of the donor star (75). We find substantially higher masses for both the black hole 
(21.2 + 2.2Mo; Fig. S8) and the donor star (40.6747 Mo). The black hole orbital semi-major 
axis derived from this reanalysis of the optical data is 0.160 + 0.013 astronomical units (au) 
(Table 1), which equates to 73 + 8 pas at our best-fitting distance of 2.22*012 kpe. This is 
consistent with the value derived directly from our VLBA astrometry (Fig. 2, Table S3). 

The higher donor mass and greater luminosity inferred from the larger distance (Table 1) 
bring the system into closer agreement with the mass-luminosity relationship for main-sequence 
hydrogen-burning stars of solar composition (7, 10). However, the measured surface composi- 
tion shows that the helium-to-hydrogen ratio is enhanced by a factor of 2.6 relative to the solar 
composition (/5). This would imply a slightly different mass-luminosity relationship if the sur- 
face abundance is indicative of the overall composition of the donor star, but remains broadly 


consistent with our revised values (Fig. 3). 

The higher mass and distance could also affect the black hole spin determined from spec- 
tral fitting of the X-ray continuum (/6). We therefore reanalysed archival X-ray data using a 
continuum-fitting method, assuming the black hole spin axis is aligned with the orbital plane. 
We find the black hole dimensionless spin parameter a, > 0.9985, close to the maximum pos- 
sible value of 1, although this could be affected by systematic uncertainties (70). Even if the 
true value is less extreme, it would still be very high, consistent with previous results derived 
from both the continuum and iron line fitting methods (/6, 77). While the spin derived from our 
analysis would be reduced if the black hole spin axis were not aligned with the orbital plane, 
even a 15° misalignment would still require a high spin, a, = 0.9696. 

As an X-ray binary system with a high-mass donor star, the black hole in Cygnus X-1 cannot 
have been spun up by accretion from its companion at the maximum theoretical rate (known as 
the Eddington limit; ~ 2 x 107" Mz yr! for Cygnus X-1), as that cannot occur for longer 
than the lifetime of the donor (~ 4 Myr for our inferred mass (/8)). The accretion time may 
be close to the age of the jet-inflated nebula surrounding the source (a few tens of kyr; (/9)). 
The current spin must therefore reflect the angular momentum of the core of the progenitor star. 
An evolutionary pathway that could explain this is main sequence mass transfer from the black 
hole progenitor to the secondary star (Case A mass transfer; (20)) with the core of the progenitor 
tidally locked and hence rapidly rotating, which can produce high black hole spins (2/). This 
evolutionary pathway for Cygnus X-1 would imply a spin axis of the black hole progenitor that 
is aligned with the orbital angular momentum. Given the very low velocity kick of 9 + 2 km s~! 
imparted to the system on black hole formation (as determined by VLBA astrometry; (22)), 
the black hole spin should still be aligned with the orbital inclination, as we assumed above. 
This scenario is also consistent with the low orbital eccentricity (6), and the lack of strong 
quasi-periodic variability seen in the X-ray power density spectra of the source (23). 

This evolutionary pathway is associated with enhanced nitrogen abundances, as observed in 
the donor spectrum (27). The transfer of enriched material from the black hole progenitor could 
also explain the high helium abundance in the spectrum of the donor star (/5). In the absence 
of convection, which is expected to be limited to a very thin surface layer in the envelope of 
a 40M, star of solar (or possibly super-solar; (/5)) metallicity, some fraction of this material 
would be retained on the surface for up to ~ 10? yr (see Supplementary Text). As a surface 
contaminant, this would not be reflective of the overall composition of the donor. Our revised 
values for the mass and luminosity of the donor star are consistent with theoretical expectations 
(Fig. 3, (70)). 

The increase in the inferred black hole mass makes Cygnus X-1 more massive than previ- 
ously observed black holes in X-ray binaries (24), surpassing M33 X-7. M33 X-7 has a substan- 
tially sub-solar metallicity (~ 0.1 times the solar metallicity, Zo; (3)). However, the metallicity 
of the mass donor in Cygnus X-1 is much higher. It has been estimated to be approximately 
twice solar (/5), although the complexity of the system makes precise measurements challeng- 
ing, and the true value may be closer to solar (70). The existence of a 21M, black hole at 
solar (or super-solar) metallicity implies that mass loss rate prescriptions (4) over-estimate the 


mass loss during the luminous blue variable or Wolf-Rayet stages of stellar evolution (25, 26). 
Assuming solar metallicity for the system, we find that either the mass loss rates in Wolf-Rayet 
winds from naked helium stars are reduced by a factor of three compared to current models (4), 
or those in luminous blue variable winds are reduced by at least a third, or both (see Supple- 
mentary Text). 

Reduced mass loss can lead to higher progenitor masses at the time of the supernova. It may 
also allow massive stars at moderate metallicities to retain hydrogen as they evolve (though 
probably not in the evolutionary history of Cygnus X-1, where the black hole’s progenitor 
likely had most of its hydrogen stripped off by the companion while the former was still on the 
main sequence; see Supplementary Text). This would change the observational signatures of 
supernovae, for example producing hydrogen-rich (pulsational) pair-instability supernovae (27). 
Enrichment of the interstellar medium by stellar winds could be reduced, and, depending on 
which phases of the stellar evolution are most impacted by reduced winds, the contribution of 
massive stars to the re-ionization of the Universe may be affected (see Supplementary Text). 

The black hole mass distribution inferred from gravitational wave events favours larger 
masses than predicted by stellar and binary evolution models. Variations in metallicity-specific 
star formation history that favour greater star formation at lower metallicity have been proposed 
as an explanation (28, 29). Reduced stellar winds would increase the mass of black holes that 
could be produced at all metallicities, such that massive gravitational wave sources could form 
at intermediate (not just very low) metallicities. This would imply that the progenitors of some 
gravitational wave events could have formed at correspondingly lower redshift, with a shorter 
delay time between binary formation and merger. 

The high spin of Cygnus X-1 (in common with most high-mass black hole X-ray binaries, 
which appear to be predominantly rapidly spinning; (30)) implies that it followed a different 
evolutionary pathway to the majority of black holes detected in gravitational wave events, which 
have spins that are either low or misaligned (3/). Given the current orbital separation, we do 
not expect Cygnus X-1 to undergo a binary black hole merger in a timescale comparable to the 
age of the Universe. 
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Figure 1: Cygnus X-1 and its best fitting model orbit on the plane of the sky. (A) astrometric 
measurements from the VLBA (red points) and the archival data from (5) (blue points). Error 
bars show the 68% confidence level. The letter labels reflect the chronological ordering of the 
observations, as detailed in Table S1. Dashed lines link the measured positions to the location 
on the fitted orbit, shown as the colored ellipse. Color bar indicates the location of the black 
hole along the line of sight, relative to the centre-of-mass of the system (shown as the black 
star), with positive values being behind the centre of mass. Arrows indicate the direction of 
orbital motion and the jet axis. (B) stacked radio image of the jet in color, with white contours 
every +( v2) times the rms noise level of 23 microJanskys (uJy) per beam. Red ellipse shows 
the synthesised beam. Although the measured positions scatter along the jet axis, the motion 
perpendicular to the jet axis is reproduced by the astrometric model (see Fig. 2). Coordinates 
are given in right ascension (RA) and declination (Dec), J2000 equinox. 
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Figure 2: Orbital displacements relative to the best-fitting one-dimensional astrometric 
model. Parallax and proper motion signatures have been subtracted. (A) The measured dis- 
placements perpendicular to the jet axis. Red points are our VLBA data, and blue points are 
the archival observations (5), with error bars showing the 6846 confidence level. Labels reflect 
the chronological ordering of the observations, as listed in Table S1. Black and magenta lines 
show 500 random draws from the posterior probability distribution of the orbital parameters, for 
radio and optical models, respectively (with the posterior medians indicated by thicker lines), 
which are consistent within the uncertainties. The data were only fitted perpendicular to the 
jet axis. (B) The measured displacements parallel to the jet axis show that the measured core 
positions are primarily downstream of the model predictions when the black hole is close to 
superior conjunction (behind the donor star; phases close to 0.0), and upstream when the black 
hole is close to inferior conjunction (phases close to 0.5), as expected for wind absorption. 
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Figure 3: Predicted mass-luminosity relations for high-mass main sequence stars. Masses 
are given in solar masses, and luminosities relative to the solar luminosity Lo. The black solid 
line shows the predicted relation for a standard composition (hydrogen mass fraction X = 0.70; 
mass fraction of heavy elements Z = 0.02 (70)). The grey dot-dashed line is for an enhanced 
helium abundance, X — 0.52, Z — 0.02 (as inferred for the surface abundance of the donor 
star; (/5)). The red solid line shows the effect of an increased metallicity, with Z = 0.03, 
X = 0.52. The magenta dotted and dashed lines show the effect of the uncertainty on the 
helium abundance (/5). The mass and luminosity determined from previous observations (6) 
are shown as the cyan triangle. The values from this work are shown as the blue circle, which 
lies closer to the theoretical relations, irrespective of composition or metallicity. Error bars 
show 68% confidence levels. 
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Parameter Median Mode Lower bound Upper bound 


i (deg) 27.51 2733 26.94 28.28 
e 0.0189 0.0186 0.0163 0.0217 
w (deg) 306.6 306.3 300.3 313.1 
M; (Mo) 406 39.8 33.5 48.3 

Ji 0.960 0.999 0.930 0.988 
TeK 31,138 31,158 30,398 31,840 
K (km s-!) 7521 75.18 74.80 75.63 
$ 0.0024 0.0023 0.0013 0.0034 
Oe 1.05 104 0.95 1.16 

Men (Mo) 212 214 18.9 234 

R, (Ro) 0*9. 2999 20.6 24.1 

log L/ Lo 5.625 — 5.606 5.547 5.698 
log(gi/cms-2) 3.348 — 3.351 3.335 3.360 
a (au) 0.244 0.243 0.231 0.256 
a, (au) 0.0838 0.0840 0.0816 0.0856 
apn (au) 0.160 0.159 0.147 0.173 


Table 1: Fitted and derived physical parameters for Cygnus X-1. The inclination 7, eccen- 
tricity e, and argument of periastron w of the orbit; the mass Mj, Roche-lobe filling factor fi, 
effective temperature 7;g, and semi-amplitude of the radial velocity curve K; for the O-star; 
a phase shift parameter @ to account for ephemeris errors; and the ratio Qot of the rotational 
frequency of the O-star to the orbital frequency (the nine parameters above the horizontal line) 
were directly fitted in the model. The black hole mass Mpy; the radius R; in solar radii Ro, 
luminosity L in solar luminosities Lo, and surface gravity gı of the O-star; and the semi-major 
axes of the full orbit a, stellar orbit a, and black hole orbit apy (the seven parameters below the 
horizontal line) were derived from the fitted parameters, using the known orbital period (/3). 
Lower and upper bounds encompass the 68% confidence interval. 
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1 Materials and Methods 


Our radio observations were made using the VLBA, as part of the Cygnus X-1 Hard-state Ob- 
servations of a Complete Binary Orbit in X-rays (CHOCBOX) project; a multi-wavelength 
observing campaign to monitor Cygnus X-1 around an entire binary orbit. 


1.1 Calibrator search 


The closest two known compact extragalactic radio calibrator sources (33) to Cygnus X-1 
are NVSS J195330+353759 (hereafter J1953+3537) and NVSS J1957404+333827 (hereafter 
J19574-3338), which lie 1.08? and 1.57? away, respectively (see Fig. S1). We chose J1957--3338 
as our primary phase reference calibrator. Despite being further from Cygnus X-1, it is less re- 
solved at our observing frequency of 8.4 GHz, allowing better calibration of the longer baselines 
to stations at Mauna Kea and St Croix, and hence better astrometric precision. 

Because systematic errors in astrometry scale with the distance between the target and the 
phase reference source (34), we conducted a calibrator search to find a closer phase reference 
source and thereby reduce the systematic uncertainties. We observed with the VLBA on 2016 
April 9, targeting 44 sources from the National Radio Astronomy Observatory (NRAO) Very 
Large Array (VLA) Sky Survey (NVSS) catalog (35) within 30' of Cygnus X-1. Observing at 
1.64 GHz, the primary beam was large enough to cover the 44 targets in 4 pointings. We used 
the multi-phase centre capability of the DiFX software correlator (36) to provide individual 
data sets for each target. Our observing bandwidth was 64 MHz, and we achieved roughly 
8 min of observations on each target source. Due to the scattering along the line of sight to the 
Cygnus region we did not use the antennas at Mauna Kea or St Croix. The brightest source 
detected (NVSS J1957544-353513, hereafter J1957+3535), at 7. 7c, was located at co-ordinates 
(32000) 19557 54:105, 4-35135'1302, 0.4? from Cygnus X-1. It appeared unresolved to our 
observations, with a peak flux density of 2.8 mJy beam" !. We therefore selected this source as 
a secondary phase reference source for our observing campaign. 


1.2 VLBA observations 


We observed Cygnus X-1 daily from 2016 May 29 through June 3, using the VLBA under 
project code BM429. Our observations covered an entire 5.6-day binary orbit, allowing us 
to track the orbital motion on the plane of the sky. Each observation was taken at a central 
frequency of 8.416 GHz, with a bandwidth of 256 MHz. We cycled between Cygnus X-1, 
our primary phase reference calibrator J1957+3338 and our secondary phase reference source 
J19574-3535, spending 100, 60 and 60s per cycle, respectively, on each source. 

Our observing runs were 12 hr in duration, with half-hour geodetic blocks at the beginning, 
middle and end of each observation. In each of these blocks, we observed a range of bright 
calibrator sources spread across the sky to solve for unmodeled tropospheric delay and clock 
errors, to improve our astrometric precision. 
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Our data were correlated using the DiFX software correlator (36), and reduced according 
to standard procedures within the 31DECI7 version of the Astronomical Image Processing 
System (AIPS; (37)). Following initial fringe fitting on J1957+3338, we stacked the data from 
all six epochs to provide a global model that was used for the final fringe fitting. The final 
phase, delay and rate solutions were interpolated to both the target source and the secondary 
phase reference calibrator, J1957+3535. To improve our amplitude calibration, we also self- 
calibrated the data on J1957+3338 from each epoch in amplitude and phase, on a 10-minute 
timescale, using our global model. We interpolated the resulting amplitude self-calibration 
solutions only (having zeroed the phases to prevent them from overriding the shorter-timescale 
phase solutions derived from the fringe fitting) to both the target source and J1957+3535. We 
then imaged both Cygnus X-1 and J1957--3535 at each epoch, and measured their positions by 
fitting the core regions with elliptical Gaussians in the image plane. Stacked images of the two 
sources made with all six epochs of data are shown in Figs. 1 and S2, respectively. 


1.3 Archival data 


Because our own data were taken over a period of just six days, they could not constrain the 
parallax and proper motion of Cygnus X-1. However, they fully sample the 5.6-day orbit, 
so when combined with archival data from 2009-2010 (5) (VLBA program code BR141), can 
measure the parallax and refine the proper motion via the extended time baseline. These archival 
data also provide an additional five samples of the binary orbit. 


1.4 Systematic uncertainties 


With a signal-to-noise > 100, our astrometric measurements were limited by systematic rather 
than statistical uncertainties. These include contributions from the ionosphere (38) and the tro- 
posphere (34), both of which scale with angular separation between target and phase calibrator 
source. 

Ionospheric wedges can produce parallax gradients across the sky, leading to systematic 
shifts of order 50 uas (v/6.7 GHz) ? per degree of offset between target and phase reference 
source, where v is the observing frequency (38). While this effect can in principle be minimized 
by fitting an ionospheric wedge using multiple phase reference sources surrounding the target, 
none of our Cygnus X-1 data sets contained sufficient well-situated calibrators to enable such 
an approach, so we adopted the above estimate of the likely systematic uncertainty. 

For our VLBA data, we used the measured offsets of the nearby secondary phase reference 
calibrator, J1957--3535, from its mean position to correct the astrometry, combining the statisti- 
cal uncertainties from each source in quadrature. This effectively phase referenced Cygnus X-1 
to J19574-3535, with the 0.4? calibrator offset implying expected systematic uncertainties from 
the troposphere of 23 and 29 uas in RA and Dec., respectively (34). We flagged local sunrise 
at each station to prevent ionospheric wedges from affecting the astrometric measurements, 
which, together with the long (12-hour) tracks, reduces the effects of any ionospheric wedge. 
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With a calibrator offset of 0.4°, the estimated ionospheric systematic uncertainty was 13 uas, 
negligible compared to the tropospheric and combined statistical uncertainties. 

The archival observations (5) used two different phase calibrator sources, J1953+3537 and 
J19574-3338, located 1.08? and 1.57? from Cygnus X-1, respectively (see Fig. S1). We used 
their position measurements (5), taking a mean of the two measurements (i.e. the positions 
measured relative to the two different phase reference sources) at each epoch to minimise the 
effect of any ionospheric wedge, and adding the statistical uncertainties from the two different 
measurements in quadrature. To this we added in quadrature both the estimated tropospheric 
systematics (38 and 47 uas in RA and Dec., respectively; (34)), and the estimated ionospheric 
systematic uncertainties of 34 uas in both dimensions. While the dependence of ionospheric 
systematic uncertainties on observing duration and time of year is not yet understood, we regard 
these estimates as conservative. The final positions of Cygnus X-1 and their uncertainties are 
given in Table S1. 

The above discussion outlines the systematic uncertainties affecting the relative astrometry 
measurements on which our results are based. To enable comparisons with other work, we also 
estimate the uncertainties affecting our absolute astrometric positions. J1957+3338 is one of the 
sources making up the third realization of the International Celestial Reference Frame (ICRF3). 
However, all our positional measurements were derived assuming the position of J1957+3338 
to be 19^5740:5499147, 33?38/277943527, which is offset by 96uas in RA and 323 as in Dec 
from the ICRF3 positions (39). The uncertainties on the ICRF3 positions are 151uas in RA 
and 181 as in Dec. The uncertainties on the absolute positions of Cygnus X-1 (Tables S1 and 
S3) and J1957--3535 (10) will be further augmented by the systematic uncertainties due to the 
offset between phase reference calibrator and target positions, as discussed above. Adding these 
contributions in quadrature, our absolute astrometric positions are uncertain by 181 as in RA 
and 375uas in Dec. 


1.5 Astrometric model fitting 


The orbital motion of the black hole can be measured in astrometric VLBI observations (5). 
Those observations showed that the sense of the black hole orbit was clockwise on the sky, but 
only measured the size of the orbit to 2c, as 0.18 + 0.09 au. Our VLBA data, corrected for the 
nominal parallax and proper motion over the six days of the campaign, showed a clockwise orbit 
that was resolved on the plane of the sky (Fig. 1). Our measured source positions are those of 
the radio core of the jet, which is presumably launched from close to the black hole. However, 
the radio photosphere is located ~ 0.5-5 mas downstream of the black hole itself (40, 41), so 
we see the projection of the true black hole orbit on the radio photosphere, potentially modified 
by orbitally- or wind-induced bending of the jet axis (40, 42). 

To sample both the orbit and proper motion of the system, we combined our VLBA data with 
the archival observations (5), and fitted the combined data set with a model incorporating proper 
motion, parallax and orbital motion. A full astrometric solution includes the reference position 
in both co-ordinates (ap and do), the proper motion in both co-ordinates (ua cos ô and us), the 
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parallax (7), plus the seven orbital parameters (orbital period P, epoch of periastron passage 
To, eccentricity e, inclination 7, semimajor axis apy, argument of periastron w, and longitude of 
the ascending node, (2). We adopted the orbital period P from the photometric ephemeris (/3). 
Constraints on the inclination of the binary orbit, the eccentricity, the argument of periastron, 
and the masses of the two components were taken from previous work (6), which, with the or- 
bital period, can be converted to a constraint on the semimajor axis via Kepler's Third Law. For 
the known inclination angle and argument of periastron (6), periastron passage (which defines 
orbital phase zero in their formalism) occurs at phase 0.109 in an alternative ephemeris (13), 
which assumed a circular orbit and defined orbital phase zero to occur at superior conjunction 
of the black hole (43). We therefore defined our reference epoch To to be the reference date for 
superior conjunction determined by (73), augmented by this phase shift. 

To fit the motion of the system on the sky, we used a Markov chain Monte Carlo (MCMC) 
approach, using the PYMC3 package (44). We adopted the Hamiltonian Monte Carlo formalism 
(HMC; (45, 46)) with a No-U-Turn Sampler (NUTS; (47)). We used Gaussian priors on ¿ and 
w based on the optical modelling work of (6). Because optical modelling cannot determine the 
sense of the orbit as seen on the plane of the sky, optically-determined inclinations of 7 and 
(180? — i) are degenerate, with values of 90° < i < 180° corresponding to the clockwise orbit 
seen in our VLBA data. Furthermore, the two components of the binary each have their own 
argument of periastron, w, and wgg, which are separated by w, = wpy + 180°. Using the 
constraints from (6), we therefore adopted a prior of 152.94 + 0.76° for the orbital inclination, 
127.6 + 5.3? for the argument of periastron of the black hole «pg, and uniform priors on six 
of the other seven parameters (œo, Ôo, Wa COS Ô, Us, 7, and apy). For the final parameter, Q, we 
assumed that the jet was aligned with the orbital angular momentum vector. The jet position 
angle is known on the plane of the sky, so with this constraint and the inclination and argument 
of periastron from the optical modelling work, we infer the longitude of the ascending node to 
be 64°. We therefore used a Gaussian prior on (2 with this mean, and a standard deviation of 1°. 
A summary of our adopted priors is given in Table S2. 

This astrometric model fitting gave a parallax of 7 = 53528 uas and an orbital semi-major 
axis for the black hole of apy = 89 + 15 uas (see Table S3 for full results). The fit residuals 
are shown in Fig. S3, and are substantially larger than we expected from our (conservative) 
estimates of the uncertainties on the data. 


1.6 Wind absorption 


To achieve a reduced x? value of 1 for their astrometric model fitting, (5) found that they had to 
adopt error floors of 80 and 160 microarcseconds (uas) in RA and Dec., respectively, which is 
again substantially larger than expected from the calibrator throw. This suggests an extra source 
of systematic uncertainty in addition to those modelled by simulations (34). 

Cygnus X-1 has a strong stellar wind, which attenuates the radio emission at frequencies 
above 2 GHz (11). The strong free-free absorption by the wind is enhanced at superior conjunc- 
tion of the black hole, which is then seen through the maximum path length through the wind. 
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At inferior conjunction, the path length is minimized and the radio emission is least absorbed. 
Because we are observing at a frequency of 8.4 GHz, where the orbital modulation of the ra- 
dio emission is 13.8 + 2.4% (40), our data are affected by wind absorption. Close to superior 
conjunction of the black hole, this would have the effect of pushing the optical depth 7 — 1 
surface out to larger distances downstream. Changes in the location of the 7 — 1 surface due 
to intrinsic variations in the electron number density or magnetic field strength in the jet would 
lead to additional, non-orbital shifts in the measured position of the radio core. 

The jet axis in Cygnus X-1 is well-known, oriented at —26° east of north (Fig. 1), and stable 
over time (5, 48). We can therefore determine the positions measured parallel and perpendic- 
ular to the jet axis, for each epoch. Because the wind absorption and intrinsic jet variability 
should only affect the parallel co-ordinate, we can conduct astrometric fits only on the per- 
pendicular component to provide an estimate of the true astrometric and orbital parameters of 
the system (32). However, when fitting in only one dimension, the reference position (o, ôo) 
and proper motion (Ha cos 0, us) become degenerate in the two sky co-ordinates. To overcome 
this problem, we used Gaussian priors for both ag and Ha cos ô, with the means and standard 
deviations taken from the posterior probability distributions of the original model fitting. 

Using this model gives a substantially smaller parallax of 7 = 458 + 35 uas (Fig. S4), and a 
semi-major axis for the black hole orbit of agy = 58 + 20 uas. Our final best-fitting parameters 
are given in Table S3. Fig. S5 shows the one-dimensional histograms of the posterior probability 
distribution of each of the nine fitted parameters, together with two-dimensional scatter plots 
showing the covariances of the parameters. The correlations between parameters arise due to 
the use of the one-dimensional fit. There is a degeneracy in position along the jet axis, as 
reflected in the mean value of the residuals parallel to the jet axis, which we found to be 86 uas 
upstream (i.e. towards the south south-east on the sky). This is well within the peak of the 
two-dimensional posterior probability distributions of (œo, 09) shown in Fig. S5. This shift has 
been corrected in Fig. 1, and also in Fig. 2, which shows the best-fitting displacements both 
perpendicular and parallel to the jet axis. 

The reduction in the amount of information available (fitting in one co-ordinate rather than 
two) increases the uncertainties in the fitted parameters. However, this model is based on well- 
characterised properties of the system, and reduces the systematic uncertainty caused by the 
radio photosphere moving up and down the jet axis as a function of orbital phase. Figs. 1 and 
2 show that, as we expected, the measured position is scattered further downstream along the 
jet axis at superior conjunction of the black hole, and back towards the black hole at inferior 
conjunction. 

Our measured parallax of 458 + 35 puas is consistent with the zero-point corrected optical 
value of ~ 470 + 40 uas from the second data release (DR2) from the Gaia space telescope 
(8, 9), which was recently refined to 468 + 15 uas in the early version of the third data release 
(eDR3; (49, 50)). Unless the uncertainty is very small, simple inversion of a measured parallax 
to determine a source distance can introduce a non-negligible bias in the distance estimation. 
With a 7.796 parallax uncertainty, we therefore adopted the standard Bayesian formalism (57) 
to convert our measured parallax to a probability density function for the source distance, using 
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an exponentially-decreasing space density prior (/2). We find a median distance of 2.22 kpc, 
with a lo range of 2.05-2.40 kpc, and a 90% Bayesian credible interval of 1.96-2.54 kpc. 


1.7 Spin-orbit misalignment 


The proper motion of Cygnus X-1 relative to its likely parent stellar association Cygnus OB3 is 
10.7+2.7kms~! (52), indicating that it formed with little to no natal kick, and possibly without 
a supernova explosion (22). Our revised distance does not change this conclusion, because it 
still falls within the distance range of 2.0 + 0.3 kpc determined for Cyg OB3 (52). 

In the absence of an asymmetric supernova kick, we would expect the angular momentum 
vector of a rapidly-spinning black hole like Cygnus X-1 (/6, 53), and hence the jet axis, to be 
aligned with the orbital angular momentum vector, which motivates our choice of a highly- 
constrained Gaussian prior for the longitude of the ascending node, Q. From the observed drift 
velocity of Cygnus X-1 relative to Cygnus OB3, we estimate the maximum asymmetric natal 
kick that the system could have received as 10-20 kms". The maximum likely misalignment 
angle should scale as the ratio of kick velocity to the pre-supernova orbital velocity. For the 
observed mass ratio between black hole and donor star, and assuming a pre-supernova orbital 
velocity of ~ 300kms~, then unless the asymmetric kick was directly opposed to the Blaauw 
kick from mass ejection (54), we estimate a maximum misalignment angle of ~ 10°. While 
inner disk misalignments of this size have been inferred from some X-ray spectral fitting results 
(55—57), other work suggests that this conclusion is sensitive to the assumed electron density in 
the accretion disk, with no misalignment being required in the case of the higher densities that 
would be appropriate for X-ray binary disks (58). 

To quantify the effect of any misalignment, we explored the effect of increasing the standard 
deviation on the Gaussian prior on Q to 15°. While the the posterior distribution for Q increased 
to 77 + 12°, all other changes were well within the uncertainties. The parallax increased by 
less than one percent, to 461 + 35 uas. The semi-major axis of the black hole orbit increased to 
63 + 20 uas, the argument of periastron increased by a degree to 126 + 5°, and the inclination 
did not change. Even for the most extreme case of a uniform prior on Q (0—360?), when the 
posterior distribution on Q increased to 95 + 18°, there was minimal further change in the 
parallax, with a median and 68% confidence interval of 464 + 35 uas. Even then, the case of 
alignment still falls within the 9096 Bayesian confidence interval. We conclude that our distance 
determination is therefore insensitive to the prior on Q, and that the data are consistent with an 
aligned jet axis. 


1.8 Updated Dynamical Model 


The increase in the distance of Cygnus X-1 compared to the previous results (5) implies that the 
absolute magnitude of the donor star is larger, in turn increasing the inferred stellar radius. The 
radius is an input to the dynamical model (6), which used optical radial velocity and light curves, 
a measurement of the projected rotational velocity of the O-star and a measurement of the radius 
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of the O-star, to determine the physical parameters of the system. Given the revised stellar radius 
arising from our distance determination, we therefore reanalysed the optical photometric and 
velocity curves of (59) and (/4) to update the estimated system parameters. 

We retained the previous stellar rotational velocity estimate (taken from (60)). The radius 
and luminosity of the O-star can be computed as a function of temperature (6), as shown in 
Fig. S6, where for simplicity we have adopted the revised distance with a conservative and 
symmetric uncertainty, 2.22 + 0.18 kpc. The larger distance implies larger radii and higher lu- 
minosities as compared to the previous results (6). We also used optical spectra (75) to improve 
previous estimates of the effective temperature, surface gravity and helium abundance of the 
O-star, which further constrain the stellar radius and luminosity. 

To model the light curve and radial velocity curve we used the ELC code (61), choosing the 
optimizer code based on the differential evolution (DE-MCMC) algorithm (62). We adopted 
the “Model D" framework from (6), which incorporates nonsynchronous rotation of the O- 
star and an eccentric orbit. We modified Model D (6) to fit for the Roche lobe filling factor 
rather than the radius of the O-star. We define this filling factor as the ratio of the distance 
between the centre of the star and the point on its surface closest to the black hole, to the 
distance between the centre of the star and the L1 Lagrange point. Our model therefore has 9 
free parameters (i, Kı, Mi, f1, Q, €, w, Qrot, Tor). We used the same table of model atmosphere 
specific intensities as (6), and hence no parameterized limb darkening law is needed. 

Because the donor star is close to Roche lobe filling (6), the Roche geometry places strong 
constraints on the surface gravity, which is close to log gı = 3.3 (6). Stellar surface gravity and 
temperature are partly degenerate in optical spectral fitting (60), so this places constraints on the 
likely effective temperature. While there are several estimates of the effective temperature of 
the donor star in the literature (/5, 60, 63, 64), we discarded (63), because it did not include the 
effects of line blanketing, which has a substantial impact on the effective temperature. Of the 
remainder (i.e., (75, 60, 64)), all suggested the effective temperature to be in the range 30,000— 
31,000 K (where we used the green spectrum of (60), as their blue spectrum gave an inconsistent 
estimate of the stellar surface gravity). We consider a broad weighted mean of these three 
effective temperature determinations, giving Tẹ = 30, 200 + 900 K. 

We used uniform priors for the 9 fitting parameters, including a broad prior on Tig of 
27,500—36,000 K. Our likelihood function is based on the x? statistic and has two parts 


2 — a2 2 
Xtotal — Xdata + Xconstraints* (S1) 


Xĉata applies to the U, B, and V light curves and the radial velocity curve, [ (6), their equation 
(2). Five additional measured properties of Cygnus X-1 are independent of orbital phase, 
namely the radius of the O-star, Ri(T.g), which is computed from the parallax for a given 
temperature; the rotational velocity of Vj, sini = 96 + 6km s^! (60); our weighted mean for 
the O-star temperature of 30, 2004-900 K; the O-star surface gravity of log gı = 3.312:0.05 (15); 
and the lack of an X-ray eclipse. A given vector of the 9 model parameters will produce values 
for all of these observed parameters, and the second part of the likelihood function is then found 
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by 


; BH WM. (Vent 06" 
Xconstraints — 6 
OR (Tog) 


Te — 30,200\*  (1ogg — 3.31 V? 
| æ — 30, 200 | og g — 3.3 ps: (S2) 
900 0.05 


where O,2 = 10° if there is an X-ray eclipse, and zero otherwise. The uncertainty of Ri (Te) 
depends on the effective temperature, hence we used c5,(7,,) in the expression above. The 
inclusion of Teg in the expression for \2,,<traints effectively imposes a Gaussian prior for that 
parameter. 


1.9 ELC Model Fitting 


To find the uncertainties on the fitted and derived parameters, we implemented the DE-MCMC 
algorithm (32). The code was run 8 times using 40 chains each. The chains were initialized 
using slightly tweaked copies of the nearly optimal model (based on the results of (6)), and each 
run had a different initial seed in the random number generator. The individual runs were done 
on computers with different CPU speeds, with between 2,500 and 26,000 generations for each 
run. 

To determine the burn-in time (how long the chains need to disperse from the initial state), 
we visually inspected plots of the chains over the first few hundred generations, finding that it 
took between about 20 and 30 generations for the chains to initially disperse. Thereafter, the 
chains reached a steady-state by generation 150 or so. To be conservative, we therefore set the 
burn-in period to be 200 generations. 

The posterior samples for the fitted and derived parameters were made by sampling the 
chains starting at generation 201 in each run, skipping enough generations to allow the chains 
to cross most of the parameter space. The eight individual posterior samples for each parameter 
for each run were combined. 

The distributions for the fitted parameters are shown in Fig. S7, and summarized with the 
derived parameters in Table 1, which gives the posterior sample median, the mode (determined 
by using 75 bins), and the lo confidence ranges. For our adopted parameters we quote the 
sample medians. 

With the exception of the inclination ¿i and the Roche lobe filling factor fi, the posterior 
probability distributions are close to symmetric (Fig. S7). In the case of the filling factor fı, the 
distribution is peaked towards the maximum value of 1.0. The median value of the distribution 
is about 0.963 and the value for the best-fitting model is 0.993. The 9596 confidence lower limit 
is 0.917. 

We find masses of 40.6* 71 Mo and 21.2 + 2.2 Mo (Fig. S8) for the O-star and black hole, 
respectively, as compared to the values of 19.2 + 1.9 Wo and 14.8 + 1.0 Mo, respectively, 
from (6). The best-fitting radius of the O-star is 22.3 + 1.8 Ro, which when combined with its 
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temperature of Teg = 31.1--0.7 KK, gives a luminosity of log( L/ Lc) = 5.63* 50. Figs. S9 and 
S10 show the light curve and radial velocity fits using these median values. We find rotation of 
the donor star consistent with being synchronous with the orbit, unlike the previous solution (6) 
which had a donor rotating at 1.4 times the orbital frequency. 

The strong dependence of the component masses on the distance results from the nearly 
Roche-lobe filling donor, which is at a known orbital period with well constrained flux and 
temperature. The angular size of the donor is set from the combination of temperature and 
luminosity. The donor radius is thus proportional to the distance to the system. Then, since 
the donor is constrained to be nearly filling its Roche lobe, and the Roche lobe size at a given 
period is set by the star's mean density, the donor mass scales as its radius cubed, and hence the 
distance cubed. This leads to the nearly doubling of the donor mass based on the 20% increase 
in distance. 

The size of the derived black hole orbit, agg = 0.160+0.013 au, corresponds to 7323-8 pas at 
our best-fitting distance of 2.22t%18 kpc. This is consistent with the size of the black hole orbit 
derived independently from our astrometric fit (58 + 20 uas), providing confidence in the result. 
This indicates that the jet cannot undergo a high degree of bending as it propagates outward 
from the black hole to the radio photosphere, or we would have measured a substantially larger 
projected orbit when fitting our VLBA astrometric data. 


1.10 Spin Measurement via the Continuum-fitting Method 


The continuum-fitting method (65) is one of the two common techniques used to estimate the 
dimensionless spin parameter of a black hole, a. = cJ/GMy, where J is the black hole 
angular momentum, c is the speed of light, and G is the gravitational constant. This method 
relies on determining the inner radius Rin of the thin accretion disk, by fitting the continuum 
X-ray spectrum to the Novikov-Thorne thin disk model (66). Rin is assumed to be located 
at the innermost stable circular orbit (ISCO) around the black hole, Rigco. Fisco is directly 
related to the dimensionless spin parameter a,, which for a prograde equatorial orbit decreases 
from 6 GM/c to GM/c? as the spin goes from a, = 0 to 1. Fitting the inner radius of the 
accretion disk then in principle constrains the black hole spin (67). However, the spins derived 
by the continuum-fitting method are sensitive to the black hole mass, the disk inclination, and 
the source distance, and the uncertainties on these parameters tend to dominate the uncertainties 
on the derived spins (65). 

The continuum fitting method was originally applied only to spectra from the high/soft X- 
ray spectral state (66, 69), to ensure that the inner radius of the disk was located at the ISCO, 
and to avoid confusion introduced by the strongly Comptonized component that is present in 
the hard and intermediate X-ray spectral states. While Cygnus X-1 does not reach the canonical 
high/soft state, (70) have shown that the inner disk radius remains within a few percent of the 
ISCO for X-ray spectra in which < 25% of the thermal seed photons are Compton upscattered. 
We applied this criterion to select six archival spectra (76, 53), and reanalysed them to constrain 
the spin parameter. We used our estimates for the distance, black hole mass and orbital inclina- 
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tion of Cygnus X-1. We used Monte Carlo simulations to quantify the uncertainty on the spin 
introduced by the combined observational uncertainties on D, M, and i (32, 55). 

Assuming that the binary orbital angular momentum vector is aligned with the spin axis of 
the black hole (as discussed above, and as expected under the evolutionary scenario discussed in 
the Supplementary text), we find that, consistent with the previous spin results (53), a, > 0.9985 
at the 30 confidence level. While the value of a, monotonically decreases with increasing 
inclination, we find that even given a misalignment angle of 15° (equivalent to 4 = 42.5°), the 
fitted spin value would still be as high as a, — 0.9696. 

While the best-fitting spin is above the maximum value expected from classical modelling 
(71), this could be attributed to systematic effects in the model arising from factors such as a 
finite disk thickness, or the inner radius not being exactly at the ISCO. Regardless, the true spin 
is likely high, as also found by studies fitting the Fe Ko line profile (17, 55—58). A detailed 
analysis of our spin-fitting results and a discussion of the inherent systematics is presented in a 
companion paper (72). 


1.11 An Updated Mass-Luminosity Relationship 


We computed a mass-luminosity relation using the Warsaw evolutionary code (52, 64). We 
made three changes relative to the model described in (7). First, we adopted a default effective 
temperature of 30,500 K (/5), which did not lead to any appreciable modification of the mass- 
luminosity relation. 

Next, we used a different method for choosing the parameter fsw (defined in (64) as a simple 
multiplicative factor applied to the stellar wind mass loss rate prescriptions of (25), to account 
for uncertainties in our knowledge of the winds in massive stars; see (73) for an alternative 
approach). Two example relations are available (7) for fixed values of fsw equal to 2 and 5 
(fixed for each mass along the relation). We calibrate the mass-luminosity relation specifically 
for the mass donor in Cygnus X-1. We use the fact that the stellar wind mass outflow rate from 
the star is known: My = —(2.57 + 0.05) x 1079 Mo yr-! (14). While this estimate does not 
include systematic uncertainties, considerations of wind clumping for comparable late O-type 
supergiants (74) suggest that these should be below a factor of two. For each stellar mass along 
the mass-luminosity sequence, we adjusted the value of f. to ensure that when the star evolves 
to reach the effective temperature 7.g = 30,500K, the calculated mass outflow rate would 
match the measured value (74). The relation was calculated for masses in the range 22-47 Mo. 
This procedure resulted in values of fsw varying from 4.15 to 0.415 as log( M / M) varied from 
1.345 to 1.668. 

Finally, we modified the assumed chemical composition. For the hydrogen mass fraction 
and metallicity, we adopted a set of default values (75) of X = 0.52 + 0.05 (corresponding to 
[He/H] = 0.42 + 0.05) and Z = 0.03 (a super-solar value set by the level of precision available 
(15)). Given the uncertainty on the metallicity, we also computed a relation for Z = 0.02. 
The assumed donor composition was found to have the largest effect on the computed mass- 
luminosity relation, as shown in Fig. 3. However, these models do not include the effects of 
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stellar rotation or tidal distortion on temperature and luminosity. 

Adopting instead the 1c upper bound on the effective temperature found above of Tig = 
31, 100K, we found that the predicted luminosity change for a given mass was < 5%. To 
account for systematic uncertainties in the assumed wind mass loss rate, we also tested the 
effect of changing the wind mass loss rates by 50%, finding the predicted luminosity change for 
a given mass to be < 3%. Our calculations are therefore insensitive to uncertainties in the wind 
mass loss rate and effective temperature. 

As shown in Fig. 3, our revised parameters of the donor star are consistent with all three 
of the assumed compositions, within the uncertainties. Should the “classical” composition of 
X = 0.70, Z = 0.02 be correct, then the true hydrogen content would be larger than the 
measured surface value of X = 0.52, possibly suggesting enrichment of the donor star during 
mass transfer (see Supplementary text). 


2 Supplementary text 


2.1 Wind mass loss rates 


Wind-driven mass loss rates of massive stars are difficult to accurately measure empirically 
and to model theoretically (see (75, 76) for reviews). We use simplified phenomenological 
prescriptions with scaling parameters to encode existing uncertainty, following the formalism 
of (4). 

We focus on two stellar evolution phases that determine the maximum mass of black holes: 
the luminous blue variable and Wolf-Rayet phases, whose wind mass loss rates are uncertain. 
Luminous blue variable winds are conjectured to prevent the star from crossing the Humphreys- 
Davidson limit (77). In our parametrised model, these are independent of metallicity (except 
indirectly, through conditions determining their onset) and eject mass at a rate 


Mtv = fiev x 107^ Ms yr !. (S3) 


Meanwhile all stars that strip themselves of their envelope through winds beyond the main 
sequence — or, in our preferred evolutionary scenario, are stripped by the companion while still 
on the main sequence — are subject to strong Wolf-Rayet winds with a metallicity-dependent 


rate of 
l LN ZNO8 
Mwr = fwr x 107 P (z) (z) Mo yr, (S4) 
© © 


where L is the star’s luminosity (78). 

We use single stellar evolution and wind models as implemented in the binary population 
synthesis code COMPAS (79, 80) to explore the impact of Wolf-Rayet and luminous blue variable 
wind prescriptions on the final black hole mass (32). We require that black holes with a mass 
equal to that inferred for Cygnus X-1 can be produced at a metallicity of Z = 0.02 = 1.4Z 6 (po- 
tentially somewhat lower than the metallicity of Cygnus X-1, so our wind reduction is conser- 
vative (/5)). In Fig. S11 we show the possible solutions for single stars in the two-dimensional 
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parameter space of fipy and fwr. We find that the mass loss rates of either luminous blue 
variable winds, or Wolf-Rayet winds, or both, must be reduced relative to the (4) defaults (cali- 
brated to the previous, lower estimates of the Cygnus X-1 black hole mass under the assumption 
that it formed at solar metallicity), which are also used in COMPAS (28). To form black holes 
with a mass consistent with that of the black hole in Cygnus X-1 from single stars at this metal- 
licity, either the mass loss rates in luminous blue variable winds have to be reduced by a third, 
or those in Wolf-Rayet winds must be reduced by two thirds. We also investigated the effect 
of the standard main sequence winds from O-stars, finding that the effect on the final remnant 
mass of switching them off entirely was marginal (when compared to the effects of varying 
either luminous blue variable or Wolf-Rayet winds), and therefore do not consider O-star winds 
further. 

Cygnus X-1 is unlikely to contain the highest-mass black hole that can be reached at its 
metallicity, so this provides only a conservative constraint on the maximum wind mass loss rate. 
This simplified analysis assumed single stars; if the progenitor of the Cygnus X-1 black hole 
was stripped of its hydrogen envelope during the main sequence as conjectured below, Wolf- 
Rayet winds will kick in earlier, again making these constraints conservative. Regardless, even 
the moderate wind strength reduction assumed here has consequences for the mass distribution 
of compact objects (8/, 82) and supernova models. For example, (83) find that varying wind as- 
sumptions, in addition to affecting the core mass, substantially impact the compactness param- 
eter of stellar cores, which is connected with the ‘explodability’ in a supernova. Reduced winds 
would also impact the enrichment of the interstellar medium by massive stars (84), and their 
role in the re-ionization of the universe; while ionising radiation from massive stars made the 
universe transparent (85), it strongly depends on stripping by winds or binary interactions (86). 


2.2 Comparison to other systems 


Dynamical black hole mass measurements of the high-mass X-ray binaries IC10 X-1 (87, 
$8) and NGC 300 X-1 (89-97) favour higher masses than our revised measurement of the 
Cygnus X-1 black hole. However, the source of the emission lines used for those measurements 
is unknown. For example (92, 93) point out that in the case of Wolf-Rayet companions, the 
emission lines likely originate in the wind and thus indicate the wind velocity rather than the 
binary orbital velocity, which would invalidate the black hole mass measurements. IC10 X-1 
and NGC 300 X-1 reside in galaxies with sub-solar metallicities, so the evolution of their pro- 
genitors is likely to differ from that of the black hole progenitor in Cygnus X-1. 

LB-1 (LS V +22 25) in the Milky Way has been suggested to contain a 70M. black 
hole (94). However, other work has challenged the interpretation of the broad Ha line used 
to determine the black hole mass, and suggested that the mass of the unobserved companion 
(which may not even be a black hole; (95)) is closer to 5M. (95-98). 
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2.3 Formation history 


We investigate a possible formation channel for Cygnus X-1. The close separation indicates 
that the binary must have interacted, with mass transfer from the progenitor of the black hole, 
the primary, onto its companion, the secondary. Such mass transfer would typically remove 
the envelope of the progenitor, which would dominate its moment of inertia and, assuming 
reasonably efficient angular momentum transport within the star, would contain the bulk of the 
angular momentum. Therefore, removing this envelope would leave behind a slowly spinning 
black hole (99—704). 

However, observations point to a rapidly rotating black hole, aligned with the orbital angular 
momentum (70). The black hole must therefore retain its birth spin, because the mass of a black 
hole would have to be doubled in order to appreciably change spin after birth (7/), which 
is inconsistent with the short lifetime of the system (705), even Eddington-limited accretion 
over the entire 4 Myr lifetime of the donor (both assumptions that are likely to be substantially 
overestimated) would only spin up the black hole to a maximum of a, = 0.3. The observed 
separation is too large for tides to spin up the stellar core of the primary after its envelope was 
removed. Chemically homogeneous evolution (1706, 107) could be responsible for the rapid spin 
of the black hole; however, given the observed mass ratio, it appears unlikely that the primary 
evolved chemically homogeneously but the secondary did not, as evidenced by its factor of — 2 
expansion from its zero-age main sequence radius. Spin-up is possible during the supernova 
itself through tidal torquing by the companion of ejecta which eventually fall back (708), but 
this generally yields lower spins than complete fallback. 

A channel for forming a rapidly rotating black hole in a close binary was proposed by (21): 
gradual case A mass transfer from the late main-sequence primary allows it to be spun up 
through tides while removing enough hydrogen-rich material to prevent post-main sequence 
expansion and loss of angular momentum through envelope stripping. While the details are 
uncertain, this channel appears to explain the similar properties of the high-mass X-ray binaries 
Cygnus X-1, LMC X-1, and M33 X-7, including orbital separations of ~ 50Rọ, and high 
spins (30) (see (709—111) for previous evolutionary studies of these systems). 

Assuming evolution through this channel, the primary would have lost some mass through 
winds during the Wolf-Rayet phase, widening the binary. However, this mass loss cannot have 
been too great, due to the mass of the black-hole remnant, and to avoid excessive spin-down. 
Eventually, the primary collapsed into an aligned, rapidly spinning black hole, likely through 
nearly complete fallback with little to no natal kick (712), consistent with the few km s! drift 
velocity of the binary relative to its Cygnus OB3 birth association (22). The lack of substantial 
mass loss and natal kick is supported by the very low binary eccentricity (70), although the fact 
that the eccentricity is inconsistent with zero may argue for a minimal mass loss of ~ 1M 
(perhaps through neutrino emission) and imperfect tidal circularisation through tides operating 
on the expanding secondary. 

The combination of case A mass transfer and Wolf-Rayet winds is presumably responsi- 
ble for the unusual surface abundances of the secondary, including surface helium abundance 
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enhanced by more than a factor of two (/5). This is different from the model of (709), who 
assumed that the surface abundances are due to partial stripping of the secondary during a pre- 
vious Roche lobe overflow episode onto the black hole; however, in their model the secondary 
is now less massive than the black hole, which is inconsistent with the present measurements. 
Enrichment during the supernova itself may have occurred in the black hole X-ray binary GRO 
J1655-40, rather than during previous mass transfer (7/5, 114). But as discussed above we 
favour nearly complete fallback in Cygnus X-1. Moreover, if the current secondary were, in 
fact, a partially stripped evolved star, we would expect it to have a higher core mass and there- 
fore a higher luminosity than a normal main sequence star of its current mass, which is not 
supported by observations (Fig. 3). 

Our scenario relies on a small quantity of enriched material being responsible for the ob- 
served secondary surface abundance. This demands that the material accreted from the primary 
is retained in a thin surface layer on the secondary, rather than being mixed throughout the star. 
Late main-sequence stars with the mass and metallicity of the secondary should only have a 
very thin convective surface layer, preventing efficient mixing (although the impact of tides and 
rotation could enhance convection). While temperature inversion in the accreted material would 
suppress the Rayleigh-Taylor instability, the more massive helium-rich material would be sus- 
ceptible to thermohaline mixing (//6, 117), which would eventually distribute the accreted ma- 
terial through the envelope, reducing the observed surface abundance. Therefore, the observed 
enhancement points to fairly recent enrichment and collapse of the black hole, likely < 10° 
years ago. This requirement is also necessary to avoid the enriched material being blown off by 
stellar winds from the secondary. This timescale is consistent with the few x 10^ yr estimate for 
the age of the jet (79), and hence the age of the black hole if the jet switched on promptly after 
the black hole was formed, as expected for wind mass transfer given the proximity of the two 
components. 

The formation scenario and its inherent uncertainties are discussed in detail elsewhere (118). 
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Figure S1: Locations of our observed sources on the plane of the sky. We show the posi- 
tions of our calibrator sources, NVSS J195330+353759 (J1953+3537), NVSS J195754+353513 
(J1957+3535), and NVSS J195740+333827 (J1957+3338), as well as that of Cygnus X-1. Hor- 
izontal and vertical axes are Right Ascension and Declination, respectively. North is upwards, 
and east is to the left. The calibrator throw from J1957+3338 is similar to both Cygnus X-1 
and the check source, J1957+3535, constraining the systematic uncertainties affecting our mea- 
surements of Cygnus X-1. All coordinates are in the International Celestial Reference System 
(ICRS). 
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Figure S2: 8.4-GHz VLBA image of J1957+3535. The image was made from a stack of 
all six epochs of data from 2016 May and June. The source is comprised of three compo- 
nents, all slightly extended relative to the VLBA beam. Contours are at levels of +2” times 
90 wy beam", where n = 0, 1,2, ... The rms noise level in the image is 26 uJy beam !. The 
fitted locations of the central component have a standard deviation of 14 and 30 microarcseconds 
(uas) in R.A. and Dec., respectively, over our six epochs of observation. The red ellipse in the 
bottom left corner represents the synthesised beam. 
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Figure S3: Residuals from the two-dimensional astrometric fit. Our model incorporated 
proper motion, parallax, and orbital motion, fitted in both Right Ascension (A) and Declination 
(B). Residuals exceed 0.1 mas, substantially larger than we expected from our estimated sys- 
tematic and statistical uncertainties, which are reflected in the sizes of our error bars, shown at 
the 68% confidence level. 
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Figure S4: The best-fitting parallax signal from our one-dimensional model fitting. We 
fitted the data perpendicular to the jet axis on the plane of the sky. (A) The parallax signal after 
removing the reference position, proper motion and orbital signatures. (B) The residuals after 
subtracting the best-fitting parallax of 7 = 458 + 35 uas. Error bars are shown at the 68% 
confidence level. 


36 


0.0023 


Relative 
probability 


0.000513 


ôo (uas) 


Bg 0142483 
458.11 3536 


—6283.33*1733 


A 7675.15 
124.761545 


5 +0.76 
152.9440-76 


eS 
aS 
E 
— ka 
yw 
DE " 
uo 
So 4 
5N 
“~~ +0.99 
$1 64.09 5.99 
S Ls —0.9t 
N Torn 
ef 
a r 
D dj 
"ENME 
5 * 
c3 
G of 
E7 q7 19.70 
FE 57.877 19.46 
bod 
S4 1 
zs PI 
ER xw H 
$ ! 
E i 
— 4 Ford 
e © tt 
5 hi 
ENI 
S 5$ tid 
i 
T T T T T T T T T T T T T T T T T T T T T d—- 
$ £ FS F P eJ NI EN EN SJ NI NI NI ov S © S 2 3 GS GS 6 4$ o 98 Q o9 d$ d$ ® OS ç QY 
SS S ES K E D 8 8 ni S gf E SK HS 9 S or gh gr MPP $ F $ s DES R 


Ha cos Ó 


(bis ur) TER Worb (°) torb (°) Qorb (°) Qorb (pas) 


ay (pas) ôo (uas) 
Figure S5: Results of our PyMC astrometric model fitting for Cygnus X-1. We fitted the 
position (a, ôo; relative to the best-fitting position, in uas), proper motions (ua Cos 0, u5; in 
uas yr- )), parallax (7, in uas), inclination angle of the orbit (i, in degrees), argument of perias- 
tron (w, in degrees), longitude of the ascending node (Q, in degrees), and semi-major axis of the 
black hole orbit (apy, in uas). Histograms show the one-dimensional posterior probability dis- 
tributions, and contour plots show the two-dimensional posterior probability distributions from 
our nine-dimensional parameter space. This shows both the spread of results and their covari- 
ances. The dashed vertical lines in the histograms represent the lo credible intervals, and the 
contour lines in the contour plots show the 1, 2 and 3o regions (from innermost to outermost, 
respectively). The red dots in the contour plots represent posterior realizations outside the 3c 
contour. 


37 


Radius (Ro) 


log(L/L) 


N 
[99] 
N 
[(9] 
CA 
[e 
CA 
CA 
N 
CA 
CA 
CA 
EN 


Temperature (kK) 


Figure S6: Stellar parameters of Cygnus X-1 as a function of its effective temperature. 
We show both the radius (A) and luminosity (B) of the O-star. The lines joining the points are 
not fits, and are shown only to highlight the broad trends in the data. The dependence of the 
bolometric corrections (see (6)) on temperature leads to the luminosity variation. 
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Figure S7: Posterior probability distributions for all nine parameters from our optical 
model fitting. N is the number of posterior samples in a given bin, out of a total of 24,360. The 
only two non-symmetric distributions are for the orbital inclination angle ? (A), and the Roche 
Lobe filling factor for the secondary, fı (E). The inclination angle is well-constrained, and the 
secondary appears very close to Roche lobe filling. We are only able to place a 9596 confidence 
lower limit on the filling factor, of > 0.914. Blue solid line shows the sample medians, and 
pink dashed lines show the lower and upper 15.8796 points, which between them encompass 
the lo parameter ranges. These values are listed in Table 1. Orange dotted lines show the 90% 
Bayesian credible intervals. 


39 


1000 


800 


400 


200 


12 14 16 18 20 22 24 26 28 
Mau (Mo) 


Figure S8: Posterior probability distribution of the black hole mass determined from our 
optical light curve and radial velocity modelling. N is the number of posterior samples in a 
given bin, out of a total of 24,360. Blue solid line shows the sample median, and pink dashed 
lines show the lower and upper 15.87% points, which between them encompass the lo range 
for the black hole mass, of 18.9—23.4 M 5. Orange dotted lines show the 90% Bayesian credible 
interval of 17.4—24.8 Ms. 
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Figure S9: Light curve results from our ELC model fitting. (A) The U-band light curve 
and the best-fitting model. (B) The U-band residuals. (C) The B-band light curve and the 
best-fitting model. (D) The B-band residuals. (E) The V-band light curve and the best-fitting 
model. (F) The V-band residuals. Error bars are shown at the 68% confidence level. The data 
are duplicated over two orbital cycles for clarity. 
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Figure S10: Radial velocity results from our ELC model fitting. (A) The phase-folded radial 
velocity measurements (data points) and the best-fitting model (red line). (B) The radial velocity 
residuals (dots). The filled circles are the averaged residuals in 30 bins, with uncertainties 
representing the standard deviation within each bin. For context, the fitted semi-amplitude of 
the radial velocity curve is 75.2 + 0.4 kms !. Data are duplicated over two orbital cycles for 
clarity. 
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Figure S11: The inferred wind strength in Cygnus X-1, relative to standard prescriptions. 
We show the maximum remnant mass from the evolution of single stars at Z — 0.02 as calcu- 
lated using COMPAS, as a function of two parameters describing the strength of luminous blue 
variable winds, fiv, and Wolf-Rayet winds, fwr. The solid and dashed white lines indicate the 
median values and 6846 confidence uncertainty on our Cygnus X-1 mass measurements. The 
star denotes default choices calibrated to the previous Cygnus X-1 mass measurement (4). Pa- 
rameter choices all fall below and to the left of it. Triangles show the isolated effects of reducing 
the mass loss rates in Wolf-Rayet winds by a factor of three, and reducing those in luminous 
blue variable winds by a third. 
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Table S1: Measured positions of the radio core of Cygnus X-1. Our images were fitted with 
an elliptical Gaussian function in the image plane (appropriately corrected for the measured 
offsets of the secondary phase reference calibrator for epoochs F-K). Tabulated uncertainties 
include both statistical and systematic uncertainties, as described in (70). Epochs A-E were 
taken in 2009-2010 under proposal code BR141 (5), and epochs F-K are our VLBA measure- 


ments from 2016, taken under proposal code BM429. 


Epoch Modified Right Ascension Uncertainty Declination Uncertainty 
Julian Date 
(days) (us) (Has) 
A 54854.703  19558"2]:6729257 4.3 35?12'57728295 59 
B 54934.574  1955872]:6728871 4.2 35*12/57721202 59 
C 55025.326 19"58™2186727746 4.2 35° 12'5!726103 58 
D 55136.023 19%58™2186726354 4.3 35° 12'5, 123796 60 
E 55221.788 19°58™2186726232 4.2 35°12'5!721701 59 
F 57537.437  19558"2]1:6706713 2.8 35?12'57682730 48 
G 57538.434  19^58"21:6706620 2.9 35° 12'5682710 50 
H 57539.432 19'58™2186706713 2.6 35° 12'5/682521 46 
I 57540.429  19^58"21:6706746 2.5 35?12/57682501 43 
J 57541.426 | 19^58"2]1:6706764 2.6 35?12/57682557 45 
K 57542.424  19^58"2]:6706679 2.6 35*12/57682604 47 


Table S2: MCMC priors adopted for the two-dimensional astrometric model fitting. 
(Qo, ôo) are the reference positions in RA and Dec. on the reference date, MJD 56198.0. U 
denotes a uniform distribution with the given range, and M denotes a normal (Gaussian) distri- 
bution with mean ju and standard deviation c. 


Parameter Description Prior distribution Units 
Qo R.A. reference position uUu(—0.2, 0.2) arcseconds 
Ôo Dec. reference position uUu(—0.2, 0.2) arcseconds 
Ha COS Ó R.A. proper motion U(—20, 20) mas yr ! 
m Dec. proper motion U(—20, 20) mas yr ! 
T Parallax U(0.1, 0.9) mas 
i Orbital inclination N (u = 152.94, o = 0.76) degrees 
GJ Argument of periastron A'(u = 127.6, 0 = 5.3) degrees 
Q Longitude of ascending node — A (us = 64.0, o = 1.0)' degrees 
QBH BH orbital semimajor axis U(0.0, 0.5) mas 


" Relative to mean of measured positions in R.A. and Dec.. 
' Taken from optical fitting results of (6). 
i Assuming jet axis and orbital angular momentum are aligned. 
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Table S3: Astrometric and orbital parameters of Cygnus X-1 from our one-dimensional 
MCMC model fitting. We used the priors described in Table S2. (a, ĉo) are the reference 
positions in RA and Dec. on the reference date, MJD 56198.0. The stated value is the median of 
the posterior distribution, with the 1c uncertainty (in brackets, applied to the final digits) defined 
as the spread between the median and percentiles 15.9 and 84.1 of the posterior probability 
distribution. The final two columns show the 90% Bayesian credible interval (the 5th and 95th 
percentiles of the posterior distribution) for the 1-D fit perpendicular to the jet axis. 


Parameter 2-D fit 1-D fit 1-D fit 1-D fit! 
5th percentile 95th percentile 

ag (19^58)* 21:6717793(17) 216717793(17) 216717770 . 2056717817 

ĉo (35° 12^) 57705435(18)  5!705497(60) 5/705399 5705597 

p, cosó (masyr~') — —3.804(5) —3.804(5) —3.801 —3.196 

[us (mas yr!) —6.312(6) —6.283(17) —6.312 —6.254 

7 (mas) 0.535(28) 0.458(35) 0.399 0.516 

(y 153.0(8) 152.9(7) 151.7 154.2 

w (°) 118(5) 125(5) 116 133 

Q (°) 64.4(1.0) 64.1(1.0) 62.4 65.7 

apy (uas) 89(15) 58(20) 25 90 


' Two-dimensional fit to the measured Right Ascension and Declination values, as de- 
scribed in (70). 

" One-dimensional fit perpendicular to the jet axis, adopting priors for o) and Ha cos ô taken 
from the posterior distribution of the standard fit, as described in (70). 

Positions given for the reference date, MJD 56198.0. 

* Inclinations of 90—180? indicate clockwise orbits. 
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